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Abstract 

From  ranging  and  tracking  applications,  to  directed  energy  weapons,  laser  tech¬ 
nology  is  very  important  to  the  military.  An  instrumental  part  of  laser  design  and 
characterization  depends  on  knowledge  of  transition  rates  between  various  atomic 
energy  levels.  Of  specific  interest,  is  determining  the  transition  probabilities  between 
the  2Ps/2  and  2P\/2  energy  levels  of  an  alkali-metal  noble-gas  laser,  which  is  found  in  a 
DPAL.  The  probabilities  of  transition,  caused  by  collisions  between  alkali  metal  and 
noble  gas  atoms,  have  previously  been  calculated  using  a  computationally  intensive 
quantum  mechanical  model.  This  research  project  aims  to  develop  a  simplified,  less 
intensive  method.  This  is  accomplished  through  the  use  of  a  semi-classical  approach, 
where  the  colliding  atoms  are  modeled  as  having  a  classical,  straight-line  trajectory. 
The  transition  probabilities  are  then  calculated  using  time-dependent  perturbation 
theory,  where  the  coupling  between  states  is  obtained  by  expressing  the  diabatic  cou¬ 
pling  potential  as  a  function  of  time.  Numerical  solutions  to  the  time-dependent 
perturbation  equations  are  obtained  for  various  forms  of  coupling  and  are  compared 
with  an  approximate  analytic  solution. 
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A  SIMPLE  MODEL  FOR  FINE  STRUCTURE  TRANSITIONS  IN 
ALKALI-METAL  NOBLE-GAS  COLLISIONS 

I.  Introduction 

Lasers  have  been  of  great  interest  to  the  Air  Force  for  many  years,  with  appli¬ 
cations  including  ranging,  tracking,  targeting,  high  power  laser  weapons,  and  coun¬ 
termeasures.  These  lasers  can  generally  be  broken  down  into  two  broad  categories 
based  on  the  lasing  medium.  These  categories  are  solid  state  and  gas  phase.  Current 
gas  phase  lasers  tend  to  have  the  benefits  of  higher  beam  quality  and  better  thermal 
properties  (heat  dissipation),  as  compared  to  solid  state  lasers.  Historically,  gas  phase 
lasers  have  been  pumped  by  chemical  reactions,  such  as  in  the  Chemical  Oxygen  Io¬ 
dine  Laser  (COIL).  The  COIL  was  invented  by  the  Air  Force  Weapons  Laboratory 
(now  Air  Force  Research  Laboratory)  in  1977,  and  was  used  aboard  the  Airborne 
Laser  (ABL). 

Current  gas-phase  laser  research  interests  include  the  Diode-Pumped  Alkali  Laser 
(DPAL),  which  has  the  advantage  of  being  electrically  pumped  by  diode  lasers.  This 
allows  the  DPAL  to  be  much  smaller,  lighter,  and  less  complex  than  the  COIL.  A 
DPAL  primarily  consists  of  a  gas  chamber,  a  laser  diode  pump,  and  mirrors  and 
other  optics.  The  gas  within  the  chamber  consists  primarily  of  a  noble  gas,  such  as 
helium,  along  with  a  relatively  small  amount  of  an  evaporated  alkali  metal,  such  as 
rubidium.  It  is  commonly  a  three-state  laser  system,  where  arrays  of  diode  lasers 
pump  electrons  in  the  alkali  metal  atoms  from  the  ground  atomic  energy  level,  2  S' 1/2, 
to  the  2 P3/2  excited  atomic  energy  level  (see  right  side  of  Figure  1).  Through  collisions 
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with  the  noble  gas  buffer  atoms,  the  excited  electrons  in  the  alkali  metal  relax  to  the 
2Pi/2  level,  creating  a  population  inversion  between  the  2P\/2  and  the  2S\/2  levels. 

This  thesis  focuses  on  the  collisional  de-excitation  of  alkali  metal  atoms  in  a  gas- 
phase  lasing  medium.  As  the  evaporated  metal  and  noble  gas  atoms  approach  each 
other  during  a  collision,  the  energy  levels  of  the  atoms  undergo  a  complex  interaction. 
The  alkali  metal’s  atomic  energy  levels,  which  are  constant  in  the  asymptotic  limit 
(at  large  atomic  separations),  become  varying  functions  of  the  atomic  separation  (R), 
as  shown  on  the  left  side  of  Figure  1. 

A  key  requirement  for  developing  and  characterizing  laser  technology  is  knowl¬ 
edge  of  the  transition  rates  of  electrons  between  their  various  atomic  energy  levels. 
Unfortunately,  a  direct  analytic  approach  to  determine  these  rates  quickly  becomes 
too  complicated  to  solve  exactly  for  all  but  the  simplest  systems.  In  response  to 
this  challenge,  numerical  techniques  have  been  developed  to  solve  more  complicated 
systems,  however  these  methods  are  computationally  intensive.  For  this  project,  a 
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Figure  1.  Adiabatic  potential  energy  surfaces  for  RbHe. 
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semi-classical  model  that  is  intended  to  reduce  the  computational  overhead  is  de¬ 
signed.  The  system  consists  of  a  collision  between  a  noble  gas  atom  and  an  alkali 
metal  atom,  and  is  intended  to  describe  the  collisional  de-excitation  in  a  gas-phase 
laser.  The  simplified  semi-classical  model  used  approximates  the  dynamics  of  the 
noble  gas  atom  using  a  straight  line  trajectory,  and  leads  to  the  development  of  a 
multivariable  rate  integral  which  is  then  used  to  determine  fine-structure  transition 
rates  of  the  alkali  atoms.  The  integration  of  this  integral  is  primarily  performed  nu¬ 
merically,  using  an  adaptive  Romberg  integration  routine  written  in  C.  An  analytic 
approximation  of  the  transition  rate  integral  is  also  described  and  compared  to  the 
numeric  methods. 

Previous  work  on  the  fine  structure  transitions  of  alkali  atoms  as  they  collide  with 
noble  gas  atoms  includes  a  full  quantum  mechanical  calculation  of  transition  rates 
based  on  wave  packet  propagation  on  potential  energy  surfaces  [5].  Comparisons  with 
these  full  quantum  mechanical  calculations  as  well  as  comparisons  with  experimental 
results,  facilitated  by  a  compilation  of  experimental  data  [3],  will  be  used  to  check 
the  validity  of  the  simplified  model. 
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II.  Theory 


2.1  The  Hamiltonian  for  a  Two- Atom  Interaction 

The  Hamiltonian  that  governs  the  dynamics  of  the  Alkali  Metal  (M)  +  Noble  Gas 
(Ng)  system  is  given  by 


H  —  ijvrad  +  Tn ■  +  Te  +  Vnn  +  feTV  +  Ke  +  V.so 


(1) 


The  individual  terms  of  this  Hamiltonian  are: 


Nuclear  Radial  Kinetic  Energy 


Tvrad  = 


Pi 


R 


‘Z^MNg 


(2) 


where  Pr  is  the  radial  momentum  operator  conjugate  to  the  internuclear  sepa¬ 


ration  R  = 


Rm  ~  RNg 


,  and  fj,MNg  is  the  nuclear  reduced  mass. 


Nuclear  Angular  Kinetic  Energy 


Tvane  - 


2  fJ'MNg 


Rm  ~  RNg 


(3) 


where  L  is  the  angular  momentum  operator,  and  Rm  and  RNg  are  the  nuclear 
position  vectors  for  the  alkali  metal  and  noble  gas  respectively. 


Electronic  Kinetic  Energy 


TP  = 


Zm~\~%N  g  p2 


E 

i— 1 


2  mP 


(4) 


where  Z  is  the  atomic  number,  and  me  is  the  mass  of  the  electron. 
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•  Nuclear- Nuclear  Potential  Energy 


Vnn  — 


M^Ng 


4vren 


Rm  —  RNg 


(5) 


where  e  is  the  magnitude  of  the  electron  charge,  and  eo  is  the  permittivity  of 
free  space 

Electronic-Nuclear  Potential  Energy 


Zm  +  Znq  o 

V'»  =  -  E  J- 

“  47Te0 

2=1 


Rm  ~  ri 


+ 


JNg 


RNg  T  i 


(6) 


where  lower  case  r  denote  electronic  position  vectors. 


Electronic-Electronic  Potential  Energy 


Vee.  = 


Zm  +  Zno  9 

e2  1 


47re0  \  fi  —  fj\ 


(7) 


Electron  Spin  Orbit  Coupling 


Vso  =  £{r)e-s  (8) 

where  £  is  the  spin-orbit  coupling  function,  i  is  the  electron’s  orbital  angular 
momentum,  and  s  is  the  electron’s  spin  angular  momentum. 

The  Hamiltonian  can  be  simplified  by  defining  the  following  combination  of  terms, 


HMNg  =  Te  +  Vnn  +  VeN  +  Vee 


(9) 
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which  leaves  the  simplified  form  of  the  Hamiltonian  for  the  M+Ng  interaction  as 


H  -  XWrad  +  fjvang  +  Hmn9  +  VSo  (10) 

2.2  The  Diabatic  Basis 


The  diabatic  basis  is  a  mixed  representation  consisting  of  the  radial  coordinate 
and  the  angular  momentum, 

R  J  J  )  (ii) 

fl  U)  / 

where  R  is  the  internuclear  separation,  J  is  the  total  angular  momentum,  is  the 
projection  of  J  onto  the  internuclear  axis,  j  is  the  electronic  angular  momentum,  and 
u  is  the  projection  of  j  onto  the  internuclear  axis.  The  2P3/2  and  2P\/2  alkali  atom 
energy  levels  correspond  to  a  total  of  six  eigenstates  in  this  basis.  With  R  omitted 
for  clarity,  the  eigenstates  are 


J  1/2 

1/2  1/2 
J  1/2 

-1/2  -1/2 


J  3/2 
3/2  3/2 


J  3/2 

1/2  1/2 


J  3/2 
-1/2  -1/2 


J  3/2  \ 

-3/2  -3/2  / 

(12) 
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For  the  M+Ng  system  considered  here,  O  =  u  in  each  of  these  eigenstates.  This  is 
because  the  projection  of  the  orbital  nuclear  angular  momentum  onto  the  internuclear 
axis  is  zero  ( J  =  L  +  j,  Jz  =  Lz  +  jz  =  jz,  with  Lz  —  0). 

Since  there  are  six  basis  vectors,  the  matrix  representation  of  the  Hamiltonian  in 
this  basis  will  yield  a  6  x  6  matrix.  Fortunately,  this  6x6  matrix  is  approximately 
block  diagonal  with  two  3x3  blocks  [5].  These  blocks  separate  the  positive  and 
negative  values  of  fl  =  u,  and  only  one  of  these  needs  to  be  considered  (see  [5]  for 
the  full  6x6  matrix  representation).  Considering  only  the  three  basis  vectors  having 
positive  momentum  projections,  the  3x3  matrix  representations  of  the  terms  of  the 
M+Ng  Hamiltonian  in  Equation  10  are  given  by 


TA+d 


-h2 

ZfJ-MNg 


J 

3/2 

3/2 

3/2 

J 

3/2 

1/2 

1/2 

J 

1/2 

1/2 

1/2 

J  3/2  \ 

J  3/2  \ 

J  1/2 

3/2  3/2  / 

1/2  1/2  / 

1/2  1/2 

d2 

dR2 

0 

0 

0 

d2 

dR2 

0 

0 

0 

d2 

dR2 

(13) 
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J 

3/2  \ 

J 

3/2  \ 

3/2 

3/2  / 

1/2 

1/2  / 

T 


Np.ne: 


—h2 

2  HMNgR2 


J(J  + 1)-| 

■[3  (J-W+l)]1/2 


J  3/2 
3/2  3/2 
J  3/2 
1/2  1/2 
J  1/2 
1/2  1/2 


(using  the  fact  that  L  =  J  —  j  and  L2  =  J2  +  j2  —  2 J  ■  j). 


-[3(</-!)G/  +  !)]1/2 


j(j  +  i)  +  f 


J  1/2 

1/2  1/2 

0 


</(</  +  1)  + 


(14) 


J 

3/2  \ 

J 

3/2  \ 

3/2 

3/2  / 

1/2 

1/2  / 

Hma^  +  Vso  — 


J  3/2 
3/2  3/2 
J  3/2 
1/2  1/2 
J  1/2 
1/2  1/2 


n  +  § 


J  1/2 

1/2  1/2 

0 


i(2£  +  n)  +  |  -3f(s-n) 


-^(s-n) 


|(E  +  2n)-e 


(15) 


»MW 


2.3  The  Adiabatic  Basis  and  Potential  Energy  Surfaces 


The  adiabatic  basis  vectors  corresponding  to  the  2Pi/2  and  2P3/2  atomic  energy 
levels  are  given  by 


In  this  basis,  the  3x3  matrix  representation  of  the  electronic  M+Ng  Hamiltonian  is 
given  by  the  diagonal  matrix 


HmJVj  +  ^  so 


(n3/2 

(pl/2 

(ni/2 


n3/2) 

n3/2 

o 

o 


o 


^1/2 

0 


Hi/s) 

0 

0 

ni/2 


(17) 


Because  there  are  no  off-diagonal  terms  in  this  matrix,  the  adiabatic  electronic  M+Ng 
Hamiltonian  does  not  contribute  to  mixing  between  eigenstates.  The  components  of 
this  matrix  have  been  previously  calculated  using  a  many-body  numerical  approach 
[1]  to  yield  adiabatic  electronic  potential  energy  surfaces,  shown  for  RbHe  in  Figure 
2.  Note  that  these  are  the  upper  states  from  the  surfaces  shown  in  Figure  1  for  a 
DPAL. 


2.4  Adiabatic-Diabatic  Transformation  and  Diabatic  Potential  Energy 
Surfaces 

As  stated,  the  potential  energy  surfaces  have  been  calculated  in  the  adiabatic  basis 
[1],  but  the  time-dependent  perturbation  theory  method  that  will  be  used  requires 
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Figure  2.  Adiabatic  potential  energy  surfaces  for  RbHe. 

the  electronic  Hamiltonian  to  be  represented  in  the  diabatic  basis.  What  is  required 
is  a  unitary  adiabatic-diabatic  transformation,  such  that 

Id-rliabat.ic  U^lHac[iabaticfJ  (1®) 

where  hJ  is  the  unitary  adiabatic-diabatic  transformation  matrix.  This  transformation 
is  found  in  [5]  by  diagonalizing  Hdiabatic  to  find  its  eigenvectors,  then  matching  those 
eigenvectors  to  the  diagonal  elements  of  Hadiabatic-  This  gives  the  elements  of  Hadiabatic 
in  terms  of  the  diagonalized  elements  of  Hdiabatic  as 

n3/2  =  n  +  | 

Si/2  =  1  +  2  (s  +  n)  +  sj  9e  +  (s  -  n)  +  4  (£  -  n)2)  (i9) 

nV2  =  |  ^-e  +  2(s  +  n)-  ^9e  +  4e(s-n)  +  4(s-n)2^) 
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These  solutions  are  then  inverted  to  obtain  the  diabatic  electronic  Hamiltonian  (given 
by  Equation  15)  in  terms  of  the  known  adiabatic  electronic  Hamiltonian  eigenvalues. 
The  transformations  are 


£  =  2111/2  nf+2Sl/2  +  \^n\/2  +  2n1/2n3/2  -  2iil/2  -  4n1/2E1/2  +  2n3/2E1/2  +  e^/2 

n  =  |  (n1/2  +  4n3/2  +  Eiy2  —  .^/n^  +  2n!/2n3/2  —  2n^2  —  ddi^E^  +  2n3/2E!/2  + 

£  —  |  n3/2  +  2n3/2  —  Ei/2  +  y/n^  +  2n!/2n3/2  —  2n^2  —  4np2£i/2  +  2n3/2E!/2  +  e^2^ 

(20) 

With  this  information,  the  diabatic  electronic  potential  energy  surfaces  can  be 
found  from  the  diagonal  terms  of  the  diabatic  electronic  Hamiltonian  matrix  (Equa¬ 
tion  15).  The  surfaces  for  RbHe  are  shown  in  Figure  3. 


Figure  3.  Diabatic  potential  energy  surfaces  for  RbHe. 
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Figure  4.  Classical  collision  between  an  alkali  Metal  atom,  M,  and  a  noble  gas  atom, 
Ng.  (The  inset  plot  displays  the  radial  coupling  function  potential  for  rubidium  and 
helium). 

2.5  Semi-Classical  Model  for  the  Probability  of  Transition  for  a  Single 
Atom 

In  order  to  simplify  the  full  quantum  mechanical  calculations  of  the  fine  structure 
transition  rates  [5],  several  approximations  are  made  to  the  Hamiltonian  matrices 
in  Equations  13,  14,  and  15.  This  is  accomplished  through  the  introduction  of  the 
semi-classical,  two-atom  collision  model.  In  this  model,  an  alkali  metal  atom  and 
a  noble  gas  atom  are  taken  to  be  point  particles.  Their  locations  are  described  by 
a  two-dimensional  Cartesian  coordinate  system  that  is  centered  on  the  alkali  metal 
and  is  stationary  in  the  alkali  metal’s  rest  frame  (see  Figure  4).  The  noble  gas  atom 
starts  at  x  =  —  oo,  and  travels  in  a  straight  line  in  the  positive  x  direction  at  a 
constant  velocity,  v.  The  collision  is  characterized  by  an  impact  parameter  b.  Due 
to  rotational  symmetry  around  the  x-axis,  only  one  plane  needs  to  be  considered, 
leaving  the  rotation  to  be  accounted  for  through  a  later  integration,  and  reducing  the 
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problem  to  a  two-dimensional  problem.  The  distance  between  the  particles  is  given 
as  R  —  \Jb2  +  (vt)2.  As  R  changes,  the  alkali  metal’s  electronic  energy  states  are 
perturbed  in  a  quantum  mechanical  interaction  through  the  radial  coupling  function, 
V ( R ).  In  this  simplified  model,  the  classical  motion  of  the  atoms  is  unaffected  by  the 
interaction. 

To  be  feasible,  this  model  must  account  for  all  terms  in  the  Hamiltonian  (Equa¬ 
tions  13,  14,  and  15).  The  tracking  of  these  terms  as  they  relate  the  full  quantum 
approach  to  the  semi-classical  approach  is  not  mathematically  rigorous,  but  is  instead 
based  on  reasonableness  arguments.  Of  course,  a  quantum  mechanical  process  cannot 
be  fully  described  in  classical  terms,  but  the  goal  is  to  distill  the  system  down  to  its 
most  important  elements  to  make  its  implementation  more  practical  and  to  make  its 
inner  workings  more  transparent.  Ultimately,  the  justification  for  any  model  is  de¬ 
termined  by  its  ability  to  predict  experimental  results.  Following  this  reasoning,  the 
classical  model  replaces  the  nuclear  radial  kinetic  energy  term  (Equation  13)  through 
its  linear  velocity  parameter.  It  also  replaces  the  diagonal  elements  of  the  nuclear 
angular  kinetic  energy  term  (Equation  14)  through  the  combination  of  linear  veloc¬ 
ity  and  impact  parameter.  This  leaves  the  off-diagonal  terms  of  the  nuclear  angular 
Hamiltonian  matrix  and  the  electronic  Hamiltonian  matrix.  Next,  note  that  the  tran¬ 
sitions  of  interest  are  between  the  J  3^2  \  state,  which  corresponds  to  the  2P3/2 

1/2  1/2  / 

atomic  energy  level,  and  the  J  \  state,  which  corresponds  to  the  2P\/2  atomic 

1/2  1/2  / 

energy  level.  The  mixing  of  these  states  is  responsible  for  the  establishment  of  a 
population  inversion  in  the  2P\/2  energy  level.  This  means  that  the  Coriolis  coupling 
within  the  2P3/2  energy  level,  between  the  '  )  and  ‘  )  states,  is  not  of 

3/2  3/2  /  1/2  1/2  / 

interest  and  is  therefore  neglected.  This  Coriolis  coupling  is  represented  by  the  off- 
diagonal  elements  of  the  nuclear  angular  kinetic  energy  Hamiltonian  (Equation  14), 
which  can  be  ignored,  now  leaving  only  the  electronic  Hamiltonian  matrix.  Without 
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J  3/2 


state.  This  allows 


the  Coriolis  coupling,  there  is  no  other  mixing  with  the 

'  3/2  3/2 

the  first  row  and  column  of  the  Hamiltonian  to  be  disregarded,  which  leaves  only  a 
2x2  matrix  that  governs  fine  structure  transitions, 


MFs  — 


J  3/2 
1/2  1/2 
J  1/2 
1/2  1/2 


J  3/2 

1/2  1/2 

i(2E  +  n)  +  | 


J  1/2 

1/2  1/2 

(e  —  n) 


-#(s-n)  §(s  +  2n)-£ 


(21) 


The  off-diagonal  elements  of  the  fine  structure  transition  matrix,  H pg,  in  Equation 


21  are  given  by 


V(R) 


(22) 


J  3/2 
1/2  1/2 


and 


J  1/2 
1/2  1/2 


where  V ( R )  is  the  radial  coupling  function,  which  mixes  the 
states.  A  positive  V(R)  is  arbitrarily  chosen,  since  the  magnitude  of  V(R)  is  the 
quantity  of  interest.  The  radial  coupling  surface  for  RbHe  is  shown  in  Figure  5.  The 


assumption  of  a  straight  line  trajectory  effectively  treats  the  diagonal  elements  of 
H pg  as  constants  and  yields  time-dependent  off-diagonal  terms, 


hfs(()  = 


(2A/2 

(2U/2 


£3/2 


2U/2) 


V(R(t))  E1/2 


(23) 


where  the  time  dependence  enters  through  R(t)  as  determined  by  the  straight  line 
trajectories.  The  time  dependent  Hamiltonian  in  Equation  23  can  be  used  with  time- 
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Figure  5.  The  radial  function  that  couples  the  adiabatic  and  diabatic  potential  energy 
surfaces  for  RbHe. 


dependent  perturbation  theory  (Section  2.6)  to  determine  the  probability  of  a  fine 
structure  transition  during  a  collision. 

2.6  Time-Dependent  Perturbation  Theory 

Paralleling  the  development  of  time-independent  perturbation  theory  in  [4],  sup¬ 
pose  there  are  two  stationary  state  wave  functions,  Vh/2,  corresponding  to  the  2P\/2 
atomic  energy  level,  and  ^3/2,  corresponding  to  the  2P3/2  atomic  energy  level.  These 
could  be  the  two  upper  states  of  a  DPAL.  Considering  just  these  two  atomic  energy 
levels,  and  a  single  electron  in  the  2P3/2  energy  level,  the  full  wave  function  is 

vb(r,  t)  =  C1/2^H2{r)e-iE^n  +  (24) 
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where  6*i/2  and  6*3/2  are  normalization  constants.  With  the  electron  starting  in  the 
upper  state,  the  probability  of  finding  the  electron  in  each  state  is 

Cy  2f  =  0 

6*3/2!  =  1 

Now,  if  these  are  truly  stationary  states,  then  this  system  is  of  no  interest  because 
transitions  between  the  2P3/2  and  2P\/2  energy  levels  are  required  to  establish  a  pop¬ 
ulation  inversion,  and  without  the  population  inversion  there  would  be  no  gain  in 
the  laser.  So  instead  consider  that  6*1/2  and  6*3/2  are  functions  of  time.  This  cor¬ 
responds  to  adding  a  time-dependent  perturbation  to  the  original  time-independent 
Hamiltonian,  H°, 

H  =  H°  +  H'{t)  (26) 

Solving  the  Schrodinger  equation  with  this  Hamiltonian  gives 


Cy  2  =  -\H'h,e-iMCV2 

6*3/2  =  H'ue-iAtCl/2  (27) 

^  _  E3/2~El/2 


which  are  exact  solutions.  To  apply  the  time-dependent  perturbation,  start  with  the 
zeroth  order  approximation 


c[%(t)  =  0 
c$W  =  1 


(28) 


which  corresponds  to  the  electron  remaining  in  its  original  state.  Next,  applying  the 
first  order  correction  by  combining  Equations  27  and  28,  and  solving  gives 


h 


f  H\  3  [t')e-lM'dt' 


0 


2  ’  2 


(29) 
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which  means  that  the  probability  that  the  electron  changed  energy  levels  is 

*  =  -f /»!.(*')*-““'<#'  (») 

o  2’2 

2.7  Determining  Transition  Rates  and  Probabilities 

When  applied  to  the  classical  model,  Equation  30  provides  the  single-atom  tran¬ 
sition  probability,  P0,  as 

OO  2 

P0(b,v)=  l-  I  V(R(b,v,t))  e~iAtdt  (31) 

—  OO 

where  the  limits  have  been  adjusted  to  model  the  noble  gas  atom  approaching  from 
a  distance  x  =  —  oo  where  V  is  negligible,  interacting  with  the  alkali  metal,  then 
proceeding  on  to  x  =  oo  where  V  is  negligible  again.  This  equation  is  the  probability 
that  an  electron  in  a  single  alkali  metal  atom  will  transition  between  electronic  energy 
levels  during  a  collision.  This  equation  only  accounts  for  a  single  collision  possessing 
a  single  velocity  and  a  single  impact  parameter.  This  probability  will  be  integrated 
across  velocity  and  impact  parameter,  but  first  some  relationships  must  be  described. 

The  experimental  data  that  the  results  of  this  thesis  will  be  compared  against 
were  recorded  as  collision  cross  sections,  <j(T),  which  are  functions  of  temperature 
[3].  In  a  hard  sphere  collision  model,  the  collision  cross  sections  can  be  found  through 
their  relationship  with  the  number  of  collisions,  N,  in  a  cylindrical  volume, 


N  =  Volume  x  Density 
—  {via)  x  n 


(32) 


where  v  is  the  mean  velocity,  t  is  the  elapsed  time,  and  n  is  the  number  of  atoms  per 
unit  volume.  Defining  a  rate  coefficient,  K,  as  the  number  of  collisions  per  unit  of 


17 


time,  per  particle  number  density  gives 


K 


N 
t  n 

V  (7 


The  collision  cross  section  can  now  be  found  through  the  relation 


(33) 


The  total  probability  is  given  by  the  ratio 


(34) 


P  =  -2- 

aQD 

K 

VQD  V 

where  aQD  is  the  quantum  defect  collision  cross  section,  defined  as  [3] 


(35) 


°QD  =  ir((r)  +r%)2 


(36) 


where  (r)  is  the  expectation  value  of  the  electron  position,  and  tr9  is  the  effective 
radius  of  the  atom. 

To  find  the  rate  coefficient,  K,  the  probability  for  the  transition  of  a  single  atom 
is  first  multiplied  by  velocity  and  then  averaged  over  velocity  using  the  Maxwell- 
Boltzmann  distribution  for  an  ideal  gas, 


K0  =  4tt 


!± 

2  nksT 


pv2  \ 

2kBT)  PQ{b,i ;)  dv 


(37) 


where  /i  is  the  reduced  mass  of  the  M+Ng  pair,  ks  is  the  Boltzmann  constant,  and  T 
is  the  temperature.  To  account  for  all  impact  parameters,  b,  an  integration  across  b 
is  performed,  where  a  factor  of  2n  is  included  to  account  for  the  cylindrical  symmetry 
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of  the  collision  geometry 


OO 


K\ 


b  K0  db 


(38) 


The  final  rate  equation  becomes 


K  =  2n  /  b  4n 


d 


27t/cbT 


3/2 


v3e 


e{~iAt)V(R(b,v,t))dt 


dv  db  (39) 


0  0  I  — oo  I 

The  rate  equation  has  units  of  length  cubed  over  time,  and  represents  the  probability 
of  transition  per  unit  time  per  particle  number  density. 


2.8  Approximations  of  the  Radial  Coupling  Function,  V(R) 

Besides  using  the  radial  coupling  function  V (R)  from  Equation  22  directly,  two 
approximate  forms  of  the  radial  coupling  function  V (R)  were  used  in  this  thesis.  The 
quantum  mechanically  calculated  V (. R )  from  Equation  22  will  now  be  referred  to  as 
Vqm(R)-  The  first  approximation,  Vrg(R),  is  a  rectangle  function  with  a  Gaussian 
onset /offset  as  shown  in  Figure  6  and  defined  as 


Vrg(R) 


Vn 


if  R<L 


Vo  exp 


( R-l ? 

2  W2 


if  R>  L 


(40) 


by  the  three  parameters  of  peak  potential  Vo,  interaction  length  L,  and  onset  width 
W.  This  “rectangular-Gaussian”  approximation  model  was  adjusted  to  be  a  close 
match  to  the  radial  coupling  potential  Vqm(R)-  An  attempt  was  made  to  find  an 
analytic  form  of  the  rate  equation  (Section  2.10)  using  this  Vrg{R)  approximation, 
but  it  was  unsuccessful.  This  led  to  the  second  form  of  approximation,  Vncit),  which 
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Figure  6.  A  slice  of  the  rectangular  Gaussian  approximation,  Vrg{R ),  to  the  coupling 
potential,  Vqm{R)-  The  potential  and  approximation  are  shown  for  RbHe,  with  impact 
parameter  6  =  0. 


is  shown  in  Figure  7  and  is  defined  as 


Vq  exp 


0 t+r )2 

2s2 


VRG(t) 


SV0 


Vq  exp 


{t-r? 

2s'2 


if  t  <  —t 
if  — r  <t<r 
if  t  >  t 


(41) 


It  is  found  by  the  transformation  of  the  explicit  independent  variable  in  Vrc{R)  from 
R  to  t  through  the  relationship, 


R2(t)  =b2 +v2t2  =>  t  =  ±  [(i?2  —  b2)/v2] 1/2  (42) 

then  defining  r  to  be  the  points  where  R  =  L 

t  =  ±[(L2  -b2)/v2]1/2  (43) 
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Figure  7.  A  slice  of  the  rectangular  Gaussian  approximation,  Vrg{R)i  to  the  coupling 
potential,  Vqm{R)-  The  potential  and  approximation  are  shown  for  RbHe,  with  v  = 
v  and  impact  parameter  b  =  0.  Vqm{R)  is  transformed  from  spatial  coordinates  to 
temporal  coordinates  according  to  the  relationship  R  =  sjb2  +  (vt)2  =  v  t. 


This  has  the  consequence  of  limiting  r  to  values  where  b  <  L,  which  introduces  an 
unknown  approximation  by  eliminating  grazing  collisions  for  which  b  >  L.  The  other 
main  difference  between  V (R)  and  V  (t)  lies  in  the  conversion  of  W  to  s,  where  W 
and  s  scale  the  widths  of  the  Gaussian  curves.  The  factor  s  is  defined  as 

W 
s  =  — 
v 

where  the  mean  velocity,  v,  is  used  rather  than  v  in  order  to  keep  this  factor  constant 
to  permit  an  analytic  integration  of  the  rate  equation,  K .  In  temporal  coordinates 
(which  the  integral  containing  V  is  in)  W  is  effectively  a  function  of  velocity,  W(v), 
while  s  is  not.  As  a  consequence,  Vrg(R)  maintains  its  similarity  to  Vqm(R),  while 
VRcit)  becomes  compressed  or  extended  as  the  relative  velocity  of  the  atoms  varies 
(see  Figures  8  and  9). 
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t  (a.u.)  xIO4 


Figure  8.  A  slice  of  the  rectangular  Gaussian  approximation,  VnG{t)t  to  the  coupling 
potential,  Vqm{R)-  The  potential  and  approximation  are  shown  for  RbHe,  with  v  = 
and  impact  parameter  b  =  0.  Notice  that  the  scaling  is  maintained  for  everything 
except  s  =  T .  Vqm(R)  is  transformed  from  spatial  coordinates  to  temporal  coordinates 
according  to  the  relationship  R  =  yjb2  +  ( vt )2  =vt. 
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Figure  9.  A  slice  of  the  rectangular  Gaussian  approximation,  to  the  coupling 

potential,  Vqm(R )•  The  potential  and  approximation  are  shown  for  RbHe,  with  v  =  2v 
and  impact  parameter  b  =  0.  Notice  that  the  scaling  is  maintained  for  everything 
except  s  =  Vqm(R)  is  transformed  from  spatial  coordinates  to  temporal  coordinates 
according  to  the  relationship  R  =  yjb2  +  ( vt )2  =vt. 
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2.9  Adiabaticity 


Adiabaticity  can  be  defined  in  different  ways,  depending  on  the  context,  so  it 
is  important  to  be  clear  what  is  meant  by  the  term.  In  this  thesis  it  is  defined  as 
the  ratio  of  two  characteristic  time  scales,  rcoiiision  and  rspin_orbit-  As  it  applies  to  the 
collision  between  two  atoms,  adiabaticity,  a,  is  defined  here  as 

^  _  Duration  of  Collision 

Spin-Orbit  Time  Scale  (45) 

_  'collision 

'spin-orbit 

The  duration  of  a  collision  is 

^collision  =  by  (46) 

where  L  is  the  interaction  length  defined  in  Section  2.5  (with  a  factor  of  two,  because 
it  is  defined  as  a  radial  interaction  length),  and  v  is  the  mean  relative  velocity  of  the 
atoms  at  a  given  temperature.  The  spin-orbit  time  scale  is  given  by 


"^spin-orbit 


2vr 

A 


(47) 


where  A  is  defined  in  Equation  27  as  the  angular  frequency  associated  with  the  2P3/2 
to  2P\/2  energy  level  transition.  Combining  these  two  definitions  provides  the  final 
form  of  the  adiabaticity  equation 


LA 
a  =  —  — 
v  7 r 


The  integrand  from  the  equation  for  the  probability  of  a  transition  in  a  single  atom, 


Equation  31,  is  given  by 


V{R(b,v,t ))  e~iAt 


(49) 
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and  is  plotted  in  Figures  10  and  11  to  illustrate  the  relationships  between  the  two 
characteristic  time  scales  for  RbHe  at  various  factors  of  mean  relative  velocity  for 
various  temperatures. 
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Figure  10.  Adiabaticity  range  for  a  collision  between  rubidium  and  helium.  Tempera¬ 
tures  range  from  300K  to  900K.  Velocities  range  from  0.59  v  to  4.0  v  (where  v  is  the 
average  relative  thermal  velocity).  The  plot  maintains  a  consistent  scaling  across  time. 


2.10  Analytic  Approximation  of  Transition  Probability 


An  analytic  approximation  for  the  total  fine  structure  transition  probability  (Equa¬ 
tion  35)  has  previously  been  derived  [8] ,  and  the  results  of  this  derivation  are  used  in 
this  thesis.  A  short  overview  of  this  analytic  approximation  is  given  here.  Combining 
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Figure  11.  Adiabaticity  range  for  a  collision  between  rubidium  and  helium.  Temper¬ 
atures  range  from  300K  to  900K.  Velocities  range  from  0.59  v  to  4.0  v  (where  v  is 
the  average  relative  thermal  velocity).  Note  that  the  scaling  of  the  axes  varies  (the 
frequency  of  the  oscillation  is  constant). 
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the  single-collision  transition  probability,  Equation  31,  and  the  definition  of  Vjic(t) 
in  Equation  41  gives 


Po{b,v) 


v 2 
io_ 

h 2 


+ 


+ 


f  f  V(t)  e~iAtdt 

—  OO 

f  dt  exp  |  —  -  iAt 

OO 
T 

f  dt  exp{—  iAt} 

-T 

o°  r  2 

f  dt  exp  |  -  iAt 


(50) 


After  defining  the  two  variables  A  and  B , 


*  =  -1/2($i)-p(-^)  (51) 

B  =  l  —  A  erfi  (^§) 
and  integrating,  the  expression  can  be  written  as, 

414? 

Po(b,v)  =  {A2cos2(Ar)  +  ABcos(Ar)sin(Ar )  +  £>2sm2(Ar)}  (52) 

nzAz 


After  summing  over  the  impact  parameter  and  taking  the  thermal  average,  as  done  in 
Section  2.7,  various  mathematical  manipulations  and  approximations  are  applied  [8] 
to  arrive  at  an  analytic  solution  for  the  rate  coefficient,  K,  and  the  final  probability 
is  found  to  be 


P  = 


I< 

aQDi 


(327 r2n2)  ( 


YlPl 


VQDV2  J  \  fL 


1C-  (A2  +  B2)^y 


(53) 


2  A/2ABexp{-%(lf)} 
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4 PIA  IIP 

'  21/4 


where  j  is  the  adiabaticity,  using  L'  =  2 L  and  /  = 
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III.  Procedure  and  Numerical  Methods 


3.1  Numerical  Systems 

For  any  numerical  routine,  it  is  critically  important  to  consider  the  numerical 
limitations  of  the  computational  system  being  used.  To  begin  the  discussion,  some 
terminology  must  be  defined.  Floating  point  numbers  are  represented  by  the  combi¬ 
nation  of  the  components, 

S  x  M  x  Be~E  (54) 

where  S  is  the  sign  bit,  M  is  the  mantissa,  B  is  the  base,  e  is  the  exponent,  and  E  is 
the  bias.  Being  a  finite  system,  each  of  these  components  has  a  set  number  of  digits 
allocated  to  its  representation.  As  a  simple,  concrete  example  to  help  visualize  some 
of  the  consequences  of  this  representation,  assume  that  we  have  a  base  10  system  with 
a  three  digit  mantissa  and  a  two  digit  exponent.  This  allows  the  mantissa  to  represent 
any  number  from  0  to  999,  and  allows  the  exponent  to  represent  any  number  from  0 
to  99.  Adding  a  bias  of  49  to  the  exponent  shifts  its  range  to  the  values  -49  through 
50.  A  number  represented  in  this  system  would  have  the  form  ±000  x  10±0°.  From 
this  example,  it  becomes  clear  how  the  number  of  digits  available  to  describe  each 
floating  point  number  limit  the  numerical  precision  and  the  range  of  the  system. 

The  relative  numerical  precision  of  the  system  (unsealed  by  the  exponent)  is  re¬ 
lated  to  the  machine  epsilon,  em,  which  is  defined  as  the  smallest  number  that  makes 
1  +  era  /  1  true.  It  is  easy  to  see  that  em  =  0.01  in  this  system,  because  this  value  of 
em  is  the  smallest  number  that  makes  the  following  true 

1  +  em  7^  1 

1  ±  0.01  ^  1  (55) 

)  ^  (100  x  10-°2) 


(100  x  10-°2)  +  (001  x  10“02 


It  is  important  to  notice  that  the  machine  epsilon,  which  is  determined  by  the 
mantissa,  is  not  the  same  as  the  smallest  number  that  the  system  can  represent,  which 
is  determined  by  the  exponent  and  bias.  In  the  example  presented  here,  em  =  1CT02, 
but  the  system  can  represent  numbers  as  small  as  10-49.  The  range  of  the  system  is 
determined  by  the  number  of  digits  in  the  exponent,  which  translates  to  the  number 
of  orders  of  magnitude  of  the  base  that  the  floating  point  representation  can  be  scaled 
by.  See  Table  1  for  a  list  of  common  system  values  used  by  most  modern  computers 
[2]- 


Table  1.  Standard  values  used  to  represent  floating  point  precision  and  double  precision 
numbers  in  modern  computers  (also  see  Equation  54).  *The  value  of  the  mantissa  can 
effectively  be  more  than  the  storage  space  allocated,  due  to  normalization  (see  [2]  for 
more  information). 


Float 

Double 

Base 

2 

2 

Mantissa  Number  of  Digits 

23(+l)* 

52(+l)* 

Exponent  Number  of  Digits 

8 

11 

1.19  x  1(T7 

2.22  x  10~16 

Minimum  Number 

1.18  x  10"38 

2.23  x  10-308 

Maximum  Number 

3.40  x  1038 

1.80  x  10308 

Rounding  errors  due  to  machine  epsilon  limitations  can  sometimes  be  minimized 
by  changing  the  order  that  calculations  are  done,  to  minimize  the  addition  or  sub¬ 
traction  of  numbers  with  highly  disparate  orders  of  magnitude.  Also,  to  avoid  errors 
that  are  caused  by  the  use  of  values  that  approach  the  machine’s  minimum  or  maxi¬ 
mum  representable  values,  appropriate  systems  of  units  can  be  chosen.  Here,  atomic 
units  were  selected  for  that  reason,  but  this  may  not  always  be  necessary  for  modern 
computer  systems,  as  long  as  double  precision  is  used. 

When  performing  a  very  large  number  of  calculations,  accumulation  of  machine 
error  can  begin  to  be  a  problem.  Assuming  that  the  error  is  not  biased  in  any 
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particular  direction,  machine  epsilon  will  accumulate  in  a  random  walk  pattern,  which 
puts  the  error  on  the  order  of  \/N em,  where  N  is  the  number  of  calculations  performed. 

3.2  Selection  of  Method  of  Quadrature 

Ideally,  the  number  of  calculations  required  in  a  numerical  routine  would  be  min¬ 
imized.  This  both  reduces  the  computation  time  and  reduces  the  chance  for  the 
accumulation  of  rounding  error.  Due  to  the  very  large  number  of  iterations  required 
to  solve  the  equations  in  this  thesis,  it  was  very  important  to  find  a  quickly  converging 
numerical  quadrature  (integration)  routine.  An  algorithm  using  only  the  trapezoidal 
rule  is  very  robust  but  can  require  many  computations  to  achieve  the  required  accu¬ 
racy  [7].  Due  to  this  large  number  of  computations,  the  accumulated  machine  error 
can  begin  to  be  significant.  The  upper  plots  in  Figures  12  and  13  show  the  expected 
convergence  of  a  trapezoidal  rule  routine  [7],  along  with  its  actual  convergence  for  a 
few  sample  cosine  integrals.  Figure  12  uses  single  floating  point  precision,  and  Figure 
13  uses  double  floating  point  precision.  Also  included  in  the  figures  is  the  approximate 
machine  accuracy  and  the  theoretical  accumulated  machine  error  for  floating  point 
numbers.  The  figures  also  show  the  expected  convergence  rate  of  Romberg  integration 
algorithms  [7],  which  are  described  below.  Note  that  the  error  accumulation  shown 
in  the  figures  is  for  a  single  integral.  If  performing  calculations  for  a  three-variable 
integral,  the  number  of  points  required  is  approximately  the  cube  of  the  number  of 
points  required  for  each  single-variable  integral,  assuming  similar  functions  for  each 
variable. 

Romberg  integration  routines  build  upon  the  strength  of  trapezoidal  rule  algo¬ 
rithms  but  converge  much  more  quickly.  The  Romberg  method  combines  multiple 
iterations  of  the  trapezoidal  rule  algorithm.  With  each  iteration  of  the  trapezoidal 
rule  algorithm,  the  integrand  is  calculated  at  an  increasing  number  of  points.  The 
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Iteration 


Figure  12.  Single-precision  convergence  for  numerical  integration  using  the  trapezoidal 
rule  (top)  and  the  Romberg  method  (bottom).  The  dashed  lines  show  the  actual  error 
in  a  few  sample  cosine  integrals,  while  the  solid  lines  show  theoretical  error.  The 
number  of  points  where  the  integrand  is  evaluated,  N,  is  related  to  the  number  of 
iterations,  j,  by  the  formula  N  =  2J_1  +  1.  The  black  dots  mark  the  points  where  the 
integrand  is  calculated  for  at  least  two  points  per  cycle.  Notice  that  both  routines 
crash  around  j  =  25  iterations.  This  is  likely  related  to  the  fact  that  this  is  where  the 
spacing  between  points  has  the  same  order  of  magnitude  as  em. 
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Iteration 

Figure  13.  Double-precision  convergence  for  numerical  integration  using  the  trapezoidal 
rule  (top)  and  the  Romberg  method  (bottom).  The  dashed  lines  show  the  actual  error 
in  a  few  sample  cosine  integrals,  while  the  solid  lines  show  theoretical  error.  The 
number  of  points  where  the  integrand  is  evaluated,  N,  is  related  to  the  number  of 
iterations,  j,  by  the  formula  N  =  2-J_1  +  1.  The  black  dots  mark  the  points  where  the 
integrand  is  calculated  for  at  least  two  points  per  cycle. 
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Romberg  method  uses  the  results  of  the  trapezoidal  rule  iterations  in  a  weighted 
extrapolation  to  determine  a  more  accurate  solution.  The  lower  plots  of  Figures  12 
and  13  show  the  actual  convergence  of  a  Romberg  routine.  This  Romberg  routine 
uses  the  same  cosine  integrals  as  used  by  the  trapezoidal  rule  algorithm  in  the  upper 
plots.  To  aid  comparison,  the  lower  plots  also  show  the  same  theoretical  information 
as  the  upper  plots.  Notice  that,  by  converging  faster,  less  accumulated  machine  error 
accrues,  and  the  Romberg  routine  is  able  to  achieve  a  couple  of  orders  of  magnitude  of 
higher  precision  as  compared  to  the  trapezoidal  rule,  before  machine  limits  dominate 
the  result.  This  effect  is  most  pronounced  in  the  plots  with  double  precision. 

Numerical  Recipes  in  C  [7]  includes  an  iterative,  adaptive  Romberg  integration 
routine  that  the  code  used  in  this  thesis  was  based  on.  For  the  first  iteration,  it 
calculates  the  value  of  the  integrand  at  the  limits  of  integration.  For  the  second 
iteration,  it  finds  the  midpoint  between  the  limits  of  integration,  and  calculates  the 
integrand  there.  It  continues  for  subsequent  iterations,  finding  the  midpoints  between 
each  pair  of  adjacent  points  from  the  previous  iteration,  and  calculating  the  integrand 
at  each  of  these  midpoints.  After  each  iteration  it  sums  the  results,  and  multiplies  by 
the  step  size,  combining  the  new  result  with  the  result  from  the  previous  iteration.  In 
this  way,  the  calculations  at  points  from  previous  iterations  are  not  wasted,  but  are 
instead  combined  with  calculations  from  new  iterations  so  that  the  solution  can  be 
refined  without  having  to  re-accomplish  previous  work.  The  total  number  of  points, 
N,  after  j  iterations  is 

Nj  =  2j"1  +  1  (56) 

So  far,  this  just  describes  an  iterative  trapezoidal  rule  algorithm.  The  Romberg 
method  expands  upon  and  refines  the  trapezoidal  rule  algorithm  by  applying  Richard¬ 
son  extrapolation  (also  known  as  Richardson’s  deferred  approach  to  the  limit).  Richard¬ 
son  extrapolation  combines  increasingly  precise  approximations  of  the  integral  in  a 
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way  that  cancels  low  order  error  terms  [2]  [7].  As  long  as  the  error  estimate  is  small, 
then  this  increased  order  of  the  error  estimate  translates  to  an  increase  in  convergence 
rate.  Richardson  extrapolation  is  based  on  the  formula 


Ain-1  _  p-1 

p  _  4  4  J±± 


(57) 


which  forms  the  table 

i? 


4° 

4° 

4 

r2 

J3 

7° 

J4 

4 

72 

J4 

73 

J4 

7° 

15 

4 

J2 

J5 

73 

J5 

(58) 


The  Romberg  method  is  seeded  by  the  7j  °,  which  are  the  results  of  successive  itera¬ 
tions  of  the  trapezoidal  rule  algorithm, 


p= o  _ 


2N, 


Nj 

-£ 

k= 1 


(f(xk+ i)  +  fM)  ~  I  =  /  f{x)  dx 


(59) 


with  each  increasing  value  of  i  being  the  next  iteration  of  the  Richardson  extrapola¬ 
tion.  Accuracy  of  values  in  the  table  increases  towards  the  bottom  and  to  the  right. 
When  the  difference  between  Pmax~1  and  I)max  is  smaller  than  the  convergence  error 
condition,  then  the  process  ends  and  is  given  as  the  solution  to  the  integral.  For 
example,  from  the  table  above, 


if  (  /|  —  Jg  <  Convergence  Error  Condition)  ,  then  set  I  =  1 5  (60) 
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where  the  (relative)  convergence  error  condition  is,  for  example, 


Convergence  Error  Condition  =  |/g|  x  10  6  (61) 

Due  to  the  oscillatory  nature  of  the  integrals  used  in  this  thesis,  another  very 
important  convergence  condition  was  imposed  in  addition  to  the  relative  error  re¬ 
quirement.  The  convergence  error  check  was  only  performed  after  the  number  of 
total  points  equaled  at  least  two  points  per  cycle  of  the  oscillating  integrand.  The 
requirement  for  having  at  least  two  points  of  integration  per  cycle  of  oscillation  is 
based  on  the  Nyquist  sampling  theorem  from  digital  signal  processing,  which  states 
that  a  continuous  wave  cannot  be  sampled  at  a  discrete  set  of  points  unless  the  sample 
rate  is  at  least  twice  the  rate  of  the  sampled  signal  [6].  To  visualize  this  requirement, 
the  large  dots  in  Figures  12  and  13  mark  the  points  where  this  condition  is  satisfied. 
Before  adding  this  condition  to  the  code,  convergence  tests  required  a  much  more 
strict  relative  error  condition,  which  in  turn  required  much  more  computation  time. 
Furthermore,  even  with  the  stricter  error  conditions,  convergence  was  very  erratic  and 
had  to  be  carefully  manually  verified.  With  this  “Nyquist  condition”  included,  the 
code  was  much  more  robust,  and  convergence  error  conditions  could  be  chosen  based 
on  finding  a  balance  between  speed  and  precision.  The  value  chosen  for  the  relative 
error  to  end  the  integration  loops  was  normally  set  to  10~6,  but  due  to  some  very 
high  oscillations  requiring  many  more  integration  points  and  in  turn  requiring  much 
more  computation  time,  chosen  values  ranged  from  10-4  to  10-6.  As  will  be  discussed 
in  more  detail,  there  were  many  considerations  and  optimizations  implemented  that 
took  the  computation  time  down  to  minutes  or  hours,  from  what  was  originally  ap¬ 
proximated  to  have  possibly  taken  months  on  a  standard  desktop  computer.  From 
what  has  already  been  presented,  it  can  be  seen  that  having  a  very  large  number 
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of  calculations  would  build  up  accumulated  error  that  could  possibly  dominate  the 
results,  bringing  any  conclusions  into  question. 

Notice  that  the  number  of  points  (see  Equation  56)  quickly  increases  as  the  number 
of  iterations  is  increased.  The  number  of  points  where  the  integrand  is  calculated 
begins  with  2  points  in  the  first  iteration,  reaches  17  points  after  5  iterations,  only  513 
points  by  10  iterations,  then  16,385  after  15  iterations,  and  a  cumbersome  1,048,577 
points  by  21  iterations.  Due  to  time  constraints,  the  code  used  in  this  thesis  was 
designed  to  report  the  inability  to  converge  if  convergence  was  not  achieved  by  25 
iterations,  or  16,777,217  points.  When  convergence  was  not  reached  by  25  iterations, 
the  convergence  error  condition  was  relaxed,  starting  from  a  relative  error  of  10~6. 
All  atomic  species  were  able  to  converge  with  at  least  a  10”4  relative  error  bound. 


3.3  Multivariable  Quadrature  -  Optimization  and  Numerical  Integration 
of  the  Rate  Equation 


The  multidimensional  integral  for  the  rate  is  given  by  Equation  39, 


K  =  27 r  /  b  Att 


n  \  3/2  f  (  _ 

1  '  /  uV  2k»T 


27 rksT 


o 


-  /  e^iAt)V{R{b,v,t))dt 


dv  db 


This  equation  is  challenging  to  integrate  for  three  primary  reasons;  (1)  because  the 
integrand  may  extend  across  many  periods  of  oscillation,  (2)  because  it  is  a  mul¬ 
tivariable  integral  with  three  variables  of  integration,  and  (3)  because  some  of  the 
integration  limits  extend  to  infinity.  The  techniques  described  below  were  used  in 
combination  with  fairly  standard  C-programming  optimization  techniques. 

The  problem  of  integrating  across  many  cycles  of  oscillation  was  solved  by  using 
the  adaptive  Romberg  integration  routine  along  with  a  sufficiently  large  number  of 
points  of  evaluation  of  the  integrand.  The  fact  that  the  Romberg  method  converges 
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relatively  quickly,  combined  with  the  fact  that  the  routine  is  adaptive,  meaning  that 
only  the  minimum  number  of  points  required  for  convergence  are  used,  is  very  impor¬ 
tant.  There  are  other  methods  of  numerical  quadrature  that  are  designed  specifically 
for  highly  oscillatory  functions,  and  a  short  detour  was  made  during  this  research 
to  see  if  they  could  be  feasible  here.  Unfortunately,  it  was  not  clear  if  the  methods 
found  could  handle  multivariable  integration  and/or  infinite  limits.  Due  to  time  con¬ 
straints  and  the  excessive  time  required  to  learn,  code,  and  test  the  various  alternate 
methods,  the  search  was  abandoned  before  a  definitive  answer  was  found.  Alternate 
methods  may  have  provided  faster  computation  times,  but  they  were  ultimately  un¬ 
necessary.  The  Romberg  method  is  very  robust  and  stable,  and  works  very  well  for 
the  calculations  performed  for  this  thesis,  as  long  as  certain  considerations  are  made, 
as  described  throughout  this  section. 


Splitting  Multidimensional  Integral  Into  Nested  Single- Variable  Inte¬ 
grals. 

To  numerically  integrate  a  multivariable  integral,  the  integral  must  first  be  broken 
up  into  separate,  nested  integrals.  It  is  already  clear  that  there  is  a  difficult  balance 
to  be  made  between  making  sure  that  the  integrand  is  calculated  at  enough  points 
to  get  an  accurate  and  precise  solution,  and  using  so  many  points  that  the  compu¬ 
tation  time  becomes  excessive  or  the  machine  error  begins  to  dominate.  One  way  to 
minimize  the  number  of  computations  required  is  by  carefully  choosing  the  way  that 
the  multivariable  integral,  Equation  39,  is  separated  into  single-variable  integrals. 
Equation  39  was  separated  as  follows: 


h(b,v)  = 


h 


e(~iAt)v(R(b,v,t))dt 


(62) 
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is  the  innermost  integral,  representing  the  probability  of  transition  for  a  single  velocity, 
v,  and  for  a  single  impact  parameter,  b.  The  middle  integral  can  be  taken  to  be 

/»oo 

l2(v)  =  2n  /  b[I3(b,v)\db  (63) 

Jo 

which  is  the  integration  across  all  impact  parameters  and  polar  angles.  This  leaves 

h  =  (d^r)  V2 1  |/2(")1  dv  (64) 

as  the  outermost  integral,  which  averages  the  product  of  the  velocity  and  probability 
across  the  Maxwell  Boltzmann  distribution. 

The  integrand  in  / l5  being  in  the  outermost  integral,  gets  evaluated  2-J1^1  +  l  times 
(see  Equation  56).  The  integrand  in  I2  must  be  evaluated  2-?2~1  +  1  times  for  each 
evaluation  of  I\,  or  (2-72-1  +  1)  (2J1_1  +  1).  The  number  of  evaluations  of  J3  is  then 
(2J3_1  +  1)  (2J2_1  +  1)  (2J1_1  +  1).  The  integral  across  t  contains  functions  of  each 
of  the  integration  variables  within  an  absolute  square,  so  it  must  be  the  innermost 
integral,  J3,  and  be  evaluated  Erst.  The  integrals  over  b  and  v  can  be  easily  separated, 
and  are  therefore  free  to  change  order  of  integration.  Since  the  innermost  of  these 
two  integrals  will  be  calculated  the  most  number  of  times,  it  should  require  the  least 
computation  time.  The  integration  across  v  is  much  more  complex  because  of  the 
Maxwell-Boltzmann  distribution,  which  led  to  the  assumption  that  it  would  require 
more  computation  time  to  evaluate,  and  should  therefore  be  the  outer  integral,  I\. 
This  leaves  the  integration  across  b  as  the  middle  integral,  I2. 

Simplifying  the  Integrals. 

The  next  step  in  preparing  the  integrals  for  numerical  integration  was  dividing 
integral  I3  into  trigonometric  functions,  by  applying  Euler’s  Formula,  eld  =  cos (9)  + 
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isin(9),  then  applying  some  algebraic  and  trigonometric  relations  to  get, 


cos(A  t)V(R(b,v,t))dt  +  /  sin(A  t)V(R(b,v,t))dt 


(65) 


which,  after  more  algebra,  becomes 


cos(A t)V(R(b,  v ,  t))dt 


+  K2 


sin(A t)V(R(b,  v,  t))dt 


(66) 


The  equations  can  be  further  simplified  and  optimized  by  combining  constants. 
Define: 


C3 

C\a 

Clb 


h 2 


2vr 

/  \3/2 

2  kgT 


(67) 


h(b,v)  =  C3  cos(Af)  V (R(b,  v,  t))dt 


+  C'i 


X^o  sin(Af)  V(R(b,v,t))dt 


(68) 


Note  that  if  V(R(b,v,t ))  is  even,  then  sin(Af)  V(R(b,v,t))dt  is  an  odd  func¬ 
tion  being  integrated  over  a  symmetric  interval,  so  the  integral  goes  to  zero,  simpli¬ 
fying  and  optimizing  the  rate  equation  even  further.  All  of  the  forms  of  V  used  in 


this  thesis  were  even,  so  the  sine  term  was  not  used. 


I2W)  —  C2 


b[I3(b,v)\db 


(69) 


roo 

h  =  cla  /  v3e^Cuv2)[I2(v)]dv  (70) 

Jo 

and  the  full  integral,  with  C4  =  C\a  C2  C3  becomes 
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poo  POO 

K  =  C4  v3e(-Clb  ^2)  /  b 

Jo  Jo 

3.4  Determining  the  Limits  of  Integration 

Possibly  the  most  difficult,  but  most  important,  optimization  technique  was  the 
determination  of  suitable  limits  of  integration.  It  relies  on  finding  the  values  of  the 
integration  variables,  beyond  which  no  meaningful  contribution  was  made  to  the  final 
result.  It  should  be  pointed  out  that  there  are  methods  of  dealing  with  infinite  limits 
that  involve  a  change  of  variable  such  as  i  — y  1/x,  but  these  methods  cause  the 
oscillations  of  the  integral  to  steadily  increase  towards  infinite  frequencies,  so  these 
methods  are  not  feasible  here.  There  are  three  sets  of  integration  limits  to  consider  in 
the  rate  equation  (Equation  39),  for  a  total  of  six  limits.  Four  of  these  limits  extend 
to  infinity,  but  fortunately  the  integrands  go  to  zero  fairly  rapidly. 

For  the  integration  across  time,  both  limits  go  to  infinity,  but  the  integrand  is 
scaled  by  the  radial  coupling  function,  V(R(t)),  which  goes  to  zero  as  t  — *  ±oo.  The 
non-negligible  range  of  these  limits  was  found  using  the  formula  xmax  =  Rmax  — 
L  +  nW,  where  n  is  an  integer,  and  the  impact  parameter  b  =  0  (6,  x,  R ,  L,  and  W 
are  described  in  sections  2.5  and  2.8).  Through  trial  and  error,  values  of  n  ranging 
from  5  to  7  were  found  to  give  relative  integration  errors  ranging  on  the  order  of  1CT6 
to  10-9.  The  limits  for  t  were  then  found  using  tumit  =  ±5nms-.  In  addition,  since 
the  integration  across  time  was  of  an  even  function  across  a  symmetric  interval,  the 
limits  were  adjusted  to  run  from  —tmax  to  zero,  and  the  result  was  multiplied  by  two. 

At  first  look,  the  required  limits  on  b  appear  to  be  similar  to  the  limits  on  x,  but 
this  is  not  necessarily  the  case.  The  rate  that  V (. R )  tapers  off  turns  out  to  have  a 
very  large  effect  on  the  time  integral,  and  any  discontinuities  from  ending  t  =  x/v 
too  early  have  a  significant  impact.  The  range  of  b  is  not  as  significant  because  the 


cos(At)  V(R(b,  v,  t))dt  ]  db  dv  (71) 
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integral  quickly  becomes  negligible  past  b  =  L  (normally).  The  adiabaticity  and 
relative  sizes  of  L  and  W  turn  out  to  have  a  significant  effect  on  the  required  b  range, 
so  care  must  be  taken  in  determining  the  limits,  especially  when  scaling  L  and  W 
(Section  5.3).  The  limit  brnax  was  set  by  the  formula  bmax  =  L+mW,  where  m  is  some 
integer.  The  limits  were  verified  by  inspection,  using  plots  similar  to  Figure  14  which 
show  the  probability  amplitude  (top)  and  scaled  absolute  square  of  the  probability 
amplitude  (bottom)  as  a  function  of  6,  for  a  sample  of  evenly  spaced  points  across 
v.  The  verification  by  inspection  consisted  of  making  sure  that  the  contribution  to 
the  total  probability  appeared  to  become  insignificant,  well  before  the  limits  were 
reached.  This  method  was  also  verified  with  a  subset  of  tests  where  the  limits  were 
extended  significantly,  and  the  result  was  compared  to  results  with  the  original  limits. 

The  velocity  limits  are  the  most  complicated,  and  choosing  them  correctly  is 
extremely  important  since  they  have  a  very  large  impact  on  computation  time,  clue 
to  their  relationship  to  the  number  of  oscillations  in  the  integrand.  The  t  limits  go  to 
infinity  as  v  goes  to  its  original  lower  limit  of  zero,  since  tumit  =  —  ¥mss.m  This  means 
that,  as  v  decreases,  the  number  of  oscillations,  and  in  turn  the  number  of  integration 
points  required,  increases  quickly.  The  upper  limit  on  v  is  important,  but  not  as 
critical.  These  limits  were  found  by  solving  for  the  values  of  v,  where  v  multiplied 
by  the  Maxwell-Boltzmann  distribution  dropped  to  1%  of  its  maximum  value.  These 
values  were  found  to  be 

Vmin  XmivV  .  . 

(72) 

^max  Xmax^ 

where  v  is  the  mean  velocity,  and  the  y  are  constants,  invariant  with  respect  to 
temperature  or  atomic  species.  The  fact  that  the  y  were  constants  greatly  simplified 
calculation  of  these  limits  for  different  temperatures  and  atomic  collision  partners. 
The  integration  limits  for  v  were  verified  by  inspection  similarly  to  the  b  limits, 
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see  Figure  15  for  example.  Due  to  the  strong  effect  of  v  on  the  time  integral  (or 
similarly,  the  effect  of  the  single-collision  adiabaticity) ,  the  limits  of  integration  over 
v  were  found  to  not  be  as  easily  described  with  the  Maxwell-Boltzmann  distribution 
as  originally  expected,  and  the  x  values  had  to  be  adjusted  accordingly.  Velocity 
limits  were  normally  set  to  vrnin  =  0.1b  and  vmax  =  4.0b.  In  hindsight,  dne  to 
the  intricate  interplay  of  all  of  the  variables  involved,  it  may  have  been  worthwhile  to 
take  the  time  to  implement  adaptive  integration  limits,  where  the  limits  automatically 
increased  until  no  additional  significant  change  in  the  final  solution  resulted. 

3.5  Calculations  Using  Vrg(^)  and  the  Analytic  Form  of  the  Rate  Equation 

Even  though  performing  calculations  using  the  analytic  form  of  the  rate  equa¬ 
tion  that  was  introduced  in  Section  2.10  is  much  faster  and  easier  than  performing 
Romberg  integration  on  the  multivariable  integral  of  Equation  39,  implementation  of 
the  analytic  approximation  is  not  completely  trivial.  This  is  because  A  in  Equation 
51  approaches  zero,  and  the  imaginary  error  function  in  B  approaches  infinity,  as  the 
quantity  As  increases.  Although  the  product  of  these  two  factors  is  a  finite  number, 
these  zeros  and  infinities  quickly  extend  beyond  the  limits  of  double  precision  float¬ 
ing  point  numbers.  Due  to  this  fact,  standard  numerical  software  such  as  Matlab  is 
unable  to  perform  calculations  for  all  of  the  atomic  elements  looked  at  in  this  thesis. 
Fortunately,  Mathematica  is  able  to  extend  beyond  the  limits  of  double  floating  point 
precision,  and  can  perform  the  required  calculations. 

3.6  Special  Considerations  for  Integration  When  Using  Vqm(R) 

When  using  Vqm(R)  (see  Section  2.8  for  description)  as  the  radial  coupling  func¬ 
tion  in  the  rate  integral,  special  considerations  had  to  be  made.  The  goal  was  to 
use  Vqm(R)  with  as  little  modification  as  possible,  but  the  prior  calculations  of  the 
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Figure  14.  Verification  of  Limits  and  Convergence  Across  Impact  Parameters.  RbHe  at  300K  is  shown. 
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Figure  15.  Verification  of  Limits  and  Convergence  Across  Velocities.  RbHe  at  300K  is  shown. 


atomic  potential  energy  surfaces  only  existed  for  R  >  3  [au]  (see  Figure  16).  To  over¬ 
come  this  difficulty,  the  choice  was  made  to  introduce  a  turning  point  in  the  noble  gas 
atom’s  trajectory  where  the  atom’s  kinetic  energy  was  equal  to  the  potential  energy 
of  the  adiabatic  potential  energy  surface.  That  raised  the  question  of  which  surface 
should  be  used,  since  the  atoms  are  in  a  linear  combination  of  different  states.  The 
choice  was  made  to  use  an  average  of  the  B and  energy  levels,  which  had  the 
benefit  of  also  resembling  the  ground  state  energy  level,  XE^.  Some  of  these  surfaces 
had  binding  regions  and  therefore  needed  to  be  modified  to  make  R  a  single-valued 
function  of  energy.  The  integration  code  used  a  lookup  table  to  find  the  two  nearest 
values  of  R  and  potential  energy,  then  performed  a  linear  interpolation  to  match  the 
given  value  of  kinetic  energy  to  the  turning  point. 

The  phase  of  the  probability  amplitude  in  Equation  31  was  set  to  zero  at  the 
point  (time)  of  impact,  and  there  was  no  phase  discontinuity  allowed  throughout  the 
duration  of  the  interaction.  This  allowed  the  integrand  in  Equation  31  to  be  easily 
split  into  pure  sine  and  cosine  terms  as  before,  simplifying  the  problem.  Note  that 
any  constant  phase  offset  applied  to  the  oscillatory  term  of  Equation  31  leaves  the 
equation  unchanged  after  taking  the  absolute  square. 

OO  2 

P'(b,v)  =  I  f  V(R(b,v,t))  e-^dt 

—  OO 

OO  2 

=  |  f  V(R(b,v,t ))  e~,At e~l(/> dt 

—  OO 

OO  2 

=  f  V(R(b,v,t))  e~lAtdt  (^) 

—  OO 

OO  2 

=  |  f  V(R(b,v,t))  e~iAtdt 

—  OO 

=  Po{b,v) 
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IV.  Results 


4.1  Results  Overview 


Table  2  summarizes  the  types  of  calculations  performed  in  this  thesis,  as  well  as 
the  collisions  pairs  for  which  calculations  were  completed.  The  table  also  summarizes 
which  data  is  compiled  from  other  sources. 


Table  2.  Summary  of  calculations  completed  and  data  compiled.  *Data  from  full  quan¬ 
tum  mechanical  calculations  covers  temperatures  up  to  approximately  400K.  ^Potas¬ 
sium  experimental  data  inconsistent.  ***Data  not  available. 
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4.2  Constants  and  Approximation  Parameters 

Table  3  lists  the  values  of  the  various  constants  and  approximations  for  all  M+Ng 
pairs  used  in  this  thesis.  Plots  of  the  quantum  mechanical  coupling  functions, 
Vqm(R),  along  with  the  approximate  coupling  functions,  Vrg(R) ,  for  all  M+Ng  pairs 
are  shown  in  Figure  16. 
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Figure  16.  Radial  Coupling  Functions  Vqm{R)  (dashed  lines)  and  Approximations  Vrg{R)  (solid  lines),  for  All  Atomic  Species. 


Table  3.  Values  used  for  atomic  data  constants  and  approximations.  *CsAr  collision 
cross  section  is  estimated. 


M+Ng 

Pair 

A 

[rad/a. u.  time] 

Reduced 
Mass,  /i 

Collision  Cross 
Section,  gqd 

L 

W 

V0 

KHe 

2.63E-04 

6.62E+03 

285.3 

4.7 

2.5 

0.0057 

KNe 

2.63E-04 

2.43E+04 

301.4 

5.15 

2.5 

0.0034 

KAr 

2.63E-04 

3.60E+04 

342.1 

6 

2.5 

0.0038 

RbHe 

1.08E-03 

6.97E+03 

306.4 

5.35 

2.5 

0.0043 

RbNe 

1.08E-03 

2.98E+04 

323.2 

5.7 

2.5 

0.00273 

RbAr 

1.08E-03 

4.96E+04 

364.2 

6.6 

2.5 

0.003 

CsHe 

2.52E-03 

7.08E+03 

338.2 

6.4 

2.5 

0.0026 

CsNe 

2.52E-03 

3.19E+04 

355.7 

6.6 

2.5 

0.0018 

CsAr 

2.52E-03 

5.60E+04 

396* 

7.6 

2.5 

0.0019 

4.3  Comparison  of  Numerically  Integrated  Vjic{t)  and  Analytic  Vrg(^) 


The  results  of  the  numerically  integrated  Vrg(^)  ate  compared  with  the  results  of 
calculations  using  the  analytic  integration  of  Vrg(^)  (see  Section  2.10)  in  Figures  17, 
18,  and  19.  This  comparison  helps  provide  validation  for  the  code,  as  well  as  for  the 
approximations  made  to  arrive  at  the  final  analytic  form.  The  comparison  is  a  near 
perfect  match  for  the  rubidium  and  cesium  pairs,  but  the  results  for  the  potassium 
pairs  differ  slightly  between  the  analytic  and  numeric  methods.  The  reason  for  this 
difference  is  likely  related  to  a  factor, 


_  A 2L2/i 
2ksT 


(74) 


that  was  introduced  in  [8]  to  simplify  the  analytic  calculations.  Xmax  was  used  to 
approximate  a  Meier  G  function  in  the  analytic  solution.  This  approximation  is  ex¬ 
pected  to  hold  for  values  of  approximately  Xmax  >  30,  which  holds  for  Rb  and  Cs,  but 
not  necessarily  for  K,  see  Table  4.  Secondly,  hypergeometric  functions  were  approx¬ 
imated  with  a  sin  term  multiplied  by  a  decaying  exponential.  This  approximation 
becomes  much  more  accurate  for  values  of  approximately  Xmax  >  100. 
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Table  4.  Values  of  approximation  parameter  Xmax  for  various  M+Ng  pairs  at  300K 
and  at  900K.  Data  in  this  thesis  was  calculated  for  temperatures  ranging  from  300K 
(higher  adiabaticity)  to  900K  (lower  adiabaticity). 


M+Ng  Pair 

A 

au] 

L 

au 

/i  [au] 

Xrnax  (300K) 

Xmax  (900K) 

KHe 

2.63E-04 

4.70 

6.62E+03 

5.32 

1.77 

KNe 

2.63E-04 

5.15 

2.43E+04 

23.4 

7.81 

KAr 

2.63E-04 

6.00 

3.60E+04 

47.2 

15.70 

RbHe 

1.08E-03 

5.35 

6.97E+03 

123 

41.0 

RbNe 

1.08E-03 

5.70 

2.98E+04 

596 

199 

RbAr 

1.08E-03 

6.60 

4.96E+04 

1330 

444 

CsHe 

2.52E-03 

6.40 

7.08E+03 

973 

324 

CsNe 

2.52E-03 

6.60 

3.19E+04 

4670 

1560 

CsAr 

2.52E-03 

7.60 

5.60E+04 

10800 

3620 
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Figure  17.  Ln(Probability)  vs  Adiabaticity  for  KHe,  KNe,  and  KAr.  Analytic  vs. 
numeric  integration  of  Frg(0- 
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Figure  18.  Ln(Probability)  vs  Adiabaticity  for  RbHe,  RbNe,  and  RbAr.  Analytic  vs. 
numeric  integration  of  Vrg (£)• 


Adiabaticity  [unitless] 

Figure  19.  Ln(Probability)  vs  Adiabaticity  for  CsHe,  CsNe,  and  CsAr.  Analytic  vs. 
numeric  integration  of  Frg(0- 
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4.4  Comparison  of  Numerically  Integrated  Vrg(£)  and  Numerically  Inte¬ 
grated  Vrg(R ) 

The  results  of  the  numerically  integrated  VRG{t)  are  compared  with  the  results  of 
the  numerically  integrated  Vrg{R )  in  Figures  20,  21,  and  22.  This  comparison  helps 
determine  the  effects  of  approximating  the  Gaussian  width  parameter  s  as  W/v.  It 
also  helps  determine  the  effect  of  limiting  the  impact  parameter  range  to  values  of 
b  <  L.  The  effect  of  integration  beyond  b  =  L  likely  had  little  significance,  based  on 
observations  of  plots  used  to  verify  the  chosen  integration  limits  (similar  to  Figure  14, 
discussed  in  Section  3.4).  Varying  the  width  of  the  Gaussian  does  have  a  significant 
impact,  and  using  a  sort  of  thermal  average  width  by  using  v  appears  to  be  a  good 
first  choice.  Although  using  v  does  not  fully  account  for  the  required  transformation 
to  go  from  Vrg(R)  to  Vrg(£),  the  shapes  of  the  curves  in  the  Ln(P )  vs.  Adiabaticity 
plots  seem  to  remain  intact  after  the  transformation.  The  v  scaling  can  easily  be 
adjusted  by  a  constant  scale  factor  (e.g.  1.5u,  etc.)  to  make  the  numerical  results 
match  experiment  or  to  make  Vrg(£)  match  Vrg{R )  (see  Section  5.3). 

4.5  Comparison  of  Numerically  Integrated  Rectangular-Gaussian  Vrg(R) 
and  Numerically  Integrated  Quantum  Mechanical  Vqm(R) 

The  results  of  the  numerically  integrated  Vrg(R)  are  compared  with  the  results 
of  the  numerically  integrated  Vqm(R)  in  Figures  23,  24,  and  25.  These  plots  raise 
what  is  possibly  the  biggest  question  that  stemmed  from  the  results  of  this  thesis. 
That  question  is,  why  did  the  radial  coupling  surfaces  that  are  the  most  likely  to  be 
the  most  correct  produce  the  highest  offset  error?  Unfortunately,  the  direct  effects  of 
replacing  the  radial  coupling  functions  with  a  square  Gaussian  approximation  are  not 
clear  here,  due  to  the  required  energy  conservation  condition  in  the  use  of  Vqm(R) 
(see  Section  3.6).  It  is  possible  to  impose  the  energy  conservation  condition  while  us- 
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Figure  20.  Ln(Probability)  vs  Adiabaticity  for  KHe,  KNe,  and  KAr.  Numeric  integra¬ 
tion  of  Vrg(0  vs.  numeric  integration  of  Vr g(-R). 


Figure  21.  Ln(Probability)  vs  Adiabaticity  for  RbHe,  RbNe,  and  RbAr.  Numeric 
integration  of  Vrg(^)  vs.  numeric  integration  of  Vrg{R)- 
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Adiabaticity  [unitless] 

Figure  22.  Ln(Probability)  vs  Adiabaticity  for  CsHe,  CsNe,  and  CsAr.  Numeric  inte¬ 
gration  of  Vf>c;(t)  vs.  numeric  integration  of  Vrg(R)- 


ing  Vrg(R).  This  was  not  fully  investigated  here,  but  some  investigation  was  made. 
As  shown  in  Figure  26  for  RbNe,  the  results  using  Vrg(R)  modified  with  energy 
conservation  approach  the  original  results  using  Vrg(R)  (without  energy  conserva¬ 
tion)  at  low  adiabaticities  (higher  velocities),  and  approach  the  results  using  Vqm{R ) 
(with  energy  conservation)  at  high  adiabaticities  (lower  velocities).  There  are  a  vast 
number  of  possible  forms  of  the  radial  coupling  function  V(R),  and  the  results  can 
vary  drastically  based  on  which  form  is  chosen.  (Efforts  in  this  thesis  eventually  be¬ 
came  focused  on  four  forms;  analytic  integration  of  Vrg(^);  numeric  integration  of 
Vrg(£),  numeric  integration  of  Vrg{R)i  and  numeric  integration  of  Vqm(R),  with  the 
approximations  matched  to  Vqm(R))- 
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Figure  23.  Ln(Probability)  vs  Adiabaticity  for  KHe,  KNe,  and  KAr.  Numeric  integra¬ 
tion  of  Vrg(R)  vs.  numeric  integration  of  Vqm(R)- 


Figure  24.  Ln(Probability)  vs  Adiabaticity  for  RbHe,  RbNe,  and  RbAr.  Numeric 
integration  of  Vrg(R)  vs.  numeric  integration  of  Vqm(R)- 
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Figure  25.  Ln(Probability)  vs  Adiabaticity  for  CsHe,  CsNe,  and  CsAr.  Numeric  inte¬ 
gration  of  Vrg{R)  vs.  numeric  integration  of  Vqm(-R). 
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Figure  26.  Ln(Probability)  vs  Adiabaticity  for  RbNe.  Investigating  the  effects  of  the 
energy  conservation  condition  by  modifying  the  numeric  integration  of  Vrg{R)  with 
energy  conservation.  Notice  that  the  results  using  the  modified  Vrg{R)  approach  the 
(unmodified)  Vr g(R)  results  at  low  adiabaticities  (higher  velocities),  and  approach  the 
Vqm(R)  results  at  high  adiabaticities  (lower  velocities). 
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V.  Discussion 


5.1  Comparison  with  Experimental  Data  -  Support  for  Validity  of  Po¬ 
tential  Energy  Surfaces 

As  mentioned  previously,  the  adiabatic  electronic  potential  energy  surfaces  were 
calculated  using  a  many-body  numerical  routine  by  Blank  [1],  These  surfaces  were 
then  converted  to  diabatic  surfaces  to  yield  adiabatic  electronic  potential  energy  sur¬ 
faces,  along  with  radial  coupling  functions  by  Lewis  [5].  Lewis  also  calculated  tem¬ 
perature  dependent  cross  sections  with  a  FORTRAN  numerical  routine  using  a  full 
quantum  mechanical  approach.  These  cross  sections  can  be  used  along  with  the 
quantum  defect  cross  sections  to  predict  the  total  probability  of  electronic  energy 
level  transition  during  a  collision  (see  Equation  35).  These  full-quantum  predictions 
can  be  used  to  validate  the  radial  coupling  surfaces  that  form  the  basis  of  the  semi- 
classical  model  in  this  thesis,  which  in  turn  can  provide  insight  into  the  effect  of  the 
approximations  made  in  this  thesis.  This  data  is  plotted  for  comparison  in  Figures  27, 
28,  and  29,  for  K,  Rb,  and  Cs  respectively.  These  figures  show  that  a  fully  quantum 
mechanical  approach  to  calculating  the  probabilities  comes  very  close  to  predicting 
the  experimental  results  for  the  potassium  pairs  and  for  the  rubidium  pairs,  although 
the  fit  is  not  as  good  for  cesium  pairs.  This  suggests  that  the  radial  coupling  func¬ 
tions  are  not  a  likely  source  for  the  discrepancy  between  the  semi- classical  models 
predictions  and  the  experimental  results. 

5.2  Choice  of  Energy  Level  for  Conservation  of  Energy  in  Vqm(R) 

As  mentioned  in  Section  3.6,  the  quantum  mechanical  (non-approximated)  radial 
coupling  functions  Vqm(R)  did  not  extend  to  R  —  0,  which  caused  the  introduction 
of  a  turning  point  based  on  conservation  of  energy.  This  led  to  the  requirement  to 
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Figure  27.  Ln(Probability)  vs  Adiabaticity  for  KHe,  KNe,  and  KAr.  Numeric  integra¬ 
tion  of  Vrg(R)  vs.  Experimental  data  and  Full  Quantum  Approach. 


Figure  28.  Ln(Probability)  vs  Adiabaticity  for  RbHe,  RbNe,  and  RbAr.  Numeric 
integration  of  Vrg(-R)  vs.  Experimental  data  and  Full  Quantum  Approach. 
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Figure  29.  Ln(Probability)  vs  Adiabaticity  for  CsHe,  CsNe,  and  CsAr.  Numeric  inte¬ 
gration  of  Vrg{R)  vs.  Experimental  data  and  Full  Quantum  Approach. 


choose  which  energy  potential  should  be  used.  The  decision  was  made  to  use  an 
average  of  the  BIh/2  and  II4/2  energy  levels.  This  average  also  had  the  benefit  of 
(loosely)  resembling  the  ground  state  energy  level,  XE^.  See  Figure  30  for  a  plot 
of  the  adiabatic  electronic  potential  energy  surfaces,  along  with  the  radial  coupling 
functions.  See  Figure  31  for  a  plot  showing  the  effect  of  the  choice  of  energy  level  for 
rubidium  and  three  different  noble  gas  species.  The  choice  of  energy  level  did  affect 
the  shape  of  the  plot,  but  it  did  not  change  the  offset  enough  to  imply  that  choosing 
the  wrong  surface  could  fully  account  for  the  discrepancy  between  the  experimental 
results  and  the  results  of  the  semi-classical  model. 
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Figure  30.  Energy  levels  of  all  M+Ng  combinations,  relative  to  their  asymptotic  limits.  An  average  between  the  BE^  and 
H1/2  energy  levels  was  chosen  to  determine  the  turning  point  for  rate  calculations  using  Vqm(R)- 
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Figure  31.  Resulting  Transition  Probabilities  for  Various  Choices  of  Conserved  Energy  Level.  RbHe  results  are  the  group  of 
lines  on  the  left,  RbNe  results  are  the  group  in  the  center,  and  RbAr  are  on  the  right  side  of  the  figure. 


5.3  Scaling  Vo,  L,  and  W  to  Make  the  Numerical  Results  for  Vrg(R)  Match 
Experimental  Results 

The  results  of  the  probability  calculations  turn  out  to  be  very  sensitive  to  the 
radial  coupling  function  approximation’s  fit  parameters  W,  Vo,  and  L.  Through 
trial  and  error,  it  was  found  that  applying  a  scale  factor  of  0.12  ~  1/8.3  to  Vo 
in  Vrg(R )  changed  the  predicted  electron  fine  structure  transition  probability  in  a 
collision  between  rubidium  and  neon  by  a  factor  of  approximately  0.014  «  1/69. 
This  brought  the  predicted  results  nearly  in  line  with  the  experimental  results.  The 
decrease  in  probability  was  even  across  the  calculated  adiabaticity  range  (temperature 
range  of  300K  to  900K).  This  Vo  scale  factor  was  then  applied  to  all  other  M+Ng 
combinations,  and  the  results  are  shown  in  Figures  33,  34,  and  35,  for  potassium, 
rubidium,  and  cesium,  respectively.  The  combined  results  for  all  M+Ng  pairs  are 
shown  in  Figure  36.  The  change  in  probability  was  approximately  1/69,  and  was  equal 
across  all  collision  pairs  and  all  adiabaticities.  The  scaled  radial  coupling  function  is 
shown  in  Figure  32. 

Through  trial  and  error,  it  was  also  found  that  applying  a  scale  factor  of  0.03  ~ 
1/33  to  L  in  Vrg(R)  changed  the  predicted  electron  fine  structure  transition  proba¬ 
bility  in  a  collision  between  rubidium  and  neon  by  an  average  factor  of  approximately 
0.035  ~  1/28.  This  aligned  the  central  values  (at  temperature  =  600K,  adiabaticity 
«  5.2)  of  the  predicted  results  and  the  experimental  results.  The  shift  in  proba¬ 
bility  depended  heavily  on  adiabaticity,  a ,  with  a  (300/1 )  decreasing  by  a  factor  of 
approximately  0.008  ~  1/126,  while  a (900/1 )  decreased  by  a  factor  of  approximately 
0.08  ~  1/12.  This  L  scale  factor  was  then  applied  to  all  other  M+Ng  combinations, 
and  the  results  are  shown  in  Figures  38,  39,  and  40,  for  potassium,  rubidium,  and 
cesium,  respectively.  The  combined  results  for  all  M+Ng  pairs  are  shown  in  Figure 
41.  The  effects  of  the  scaled  L  varied  greatly  across  collision  pair,  as  well  as  across 
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Figure  32.  Vrg(-R)  radial  coupling  function  with  fit  parameter,  Vq,  scaled  by  a  factor  of 
0.12.  Data  is  shown  for  RbHe,  and  Vqm(R)  is  plotted  for  comparison. 
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Figure  33.  Scaled  Vrq (_R)  radial  coupling  fit  parameter,  Vq,  for  KHe,  KNe,  and  KAr. 
Experimental  data  is  plotted  for  comparison.  Vq  is  scaled  by  a  factor  of  0.12. 


Figure  34.  Scaled  Vrg(R)  radial  coupling  fit  parameter,  Vo?  f°r  RbHe,  RbNe,  and  RbAr. 
Experimental  data  is  plotted  for  comparison.  Vq  is  scaled  by  a  factor  of  0.12. 
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Figure  35.  Scaled  Vrg(R)  radial  coupling  fit  parameter,  Vq,  for  CsHe,  CsNe,  and  CsAr. 
Experimental  data  is  plotted  for  comparison.  Vq  is  scaled  by  a  factor  of  0.12. 


adiabaticity,  with  changes  ranging  from  a  factor  of  1/2.4  for  KNe  at  300K,  to  1/187 
for  CsAr  at  300K.  The  scaled  radial  coupling  function  is  shown  in  Figure  37. 

It  was  then  also  found,  through  trial  and  error,  that  applying  a  scale  factor  of  2.6 
to  W  in  Vrg{R)  decreased  the  predicted  electron  fine  structure  transition  probability 
in  a  collision  between  rubidium  and  neon  by  an  average  factor  of  approximately 
0.013  ~  1/78.  This  brought  the  predicted  results  for  RbNe  nearly  in  line  with  the 
experimental  results.  The  shift  in  probability  depended  on  adiabaticity,  with  a/300/1 ) 
decreasing  by  a  factor  of  approximately  0.018  ~  1/55,  while  a(900/l)  decreased  by 
a  factor  of  approximately  0.009  ~  1/111.  This  scale  factor  was  then  applied  to  all 
other  M+Ng  combinations,  and  the  results  are  shown  in  Figures  43,  44,  and  45,  for 
potassium,  rubidium,  and  cesium,  respectively.  The  results  varied  across  collision 
pairs  and  adiabaticities,  but  the  variation  was  not  as  large  as  it  was  for  the  scaling 
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Figure  36.  Scaled  Vr g{R)  radial  coupling  fit  parameter,  Vo,  for  all  M+Ng  pairs.  Experimental  data  is  plotted  for  comparison. 
Vq  is  scaled  by  a  factor  of  0.12. 


Figure  37.  Vrg(R)  radial  coupling  function  with  fit  parameter,  L,  scaled  by  a  factor  of 
0.03.  Data  is  shown  for  RbHe,  and  Vqm(R)  is  plotted  for  comparison. 
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Figure  38.  Scaled  Vrg(R)  radial  coupling  fit  parameter,  L,  for  KHe,  KNe,  and  KAr. 
Experimental  data  is  plotted  for  comparison.  L  is  scaled  by  a  factor  of  0.03. 
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Figure  39.  Scaled  Vrg(R)  radial  coupling  fit  parameter,  L,  for  RbHe,  RbNe,  and  RbAr. 
Experimental  data  is  plotted  for  comparison.  L  is  scaled  by  a  factor  of  0.03. 
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Figure  40.  Scaled  Vrg{R)  radial  coupling  fit  parameter,  L,  for  CsHe,  CsNe,  and  CsAr. 
Experimental  data  is  plotted  for  comparison.  L  is  scaled  by  a  factor  of  0.03. 


of  L.  The  combined  results  for  all  M+Ng  pairs  are  shown  in  Figure  46.  The  scaled 
radial  coupling  function  is  shown  in  Figure  42. 

It  is  very  interesting  to  note  that  different  scale  factors  shift  the  transition  proba¬ 
bilities  of  specific  elements  by  different  amounts,  depending  on  the  temperature  (and 
therefore  depending  on  the  adiabaticity).  This  means  that  it  should  be  possible  to  fit 
the  predicted  transition  probabilities  across  a  range  of  temperature  (and  adiabaticity) 
values,  rather  than  just  attempting  to  minimize  an  average  error.  This  is  currently 
easier  said  than  done,  due  to  the  fact  that  the  best-fit  scale  factors  are  currently 
found  through  guessing  a  value,  performing  a  test  run,  checking  the  results,  and  then 
making  another  guess  until  the  fit  is  satisfactory.  Having  the  analytic  approximation 
of  Vrg(^)  may  make  this  task  easier. 
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Figure  41.  Scaled  Vrg(R)  radial  coupling  fit  parameter,  L,  for  all  M+Ng  pairs.  Experimental  data  is  plotted  for  comparison. 
L  is  scaled  by  a  factor  of  0.03. 


Figure  42.  Vf ig{R)  radial  coupling  function  with  fit  parameter,  W,  scaled  by  a  factor  of 
2.6.  Data  is  shown  for  RbHe,  and  Vqm(R)  is  plotted  for  comparison. 
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Figure  43.  Scaled  Vrg(R)  radial  coupling  fit  parameter,  W,  for  KHe,  KNe,  and  KAr. 
Experimental  data  is  plotted  for  comparison.  W  is  scaled  by  a  factor  of  2.6. 


Figure  44.  Scaled  Vj ig{R)  radial  coupling  fit  parameter,  W,  for  RbHe,  RbNe,  and  RbAr. 
Experimental  data  is  plotted  for  comparison.  W  is  scaled  by  a  factor  of  2.6. 
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Figure  45.  Scaled  Vrg(R)  radial  coupling  fit  parameter,  W,  for  CsHe,  CsNe,  and  CsAr. 
Experimental  data  is  plotted  for  comparison.  W  is  scaled  by  a  factor  of  2.6. 
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Figure  46.  Scaled  Vrg{R)  radial  coupling  fit  parameter,  W,  for  all  M+Ng  pairs.  Experimental  data  is  plotted  for  comparison. 
W  is  scaled  by  a  factor  of  2.6. 


VI.  Conclusion 


6.1  Conclusion 

There  are  a  few  observations  that  stand  out  from  the  results  of  the  calculations 
based  on  the  semi-classical  model  presented  in  this  thesis.  One  of  the  disappointing 
observations  is  that  the  calculations  using  the  semi-classical  model  are  off  by  an  offset 
ranging  from  about  e3  ~  20  (using  VrG(£))  to  e8  ~  3,000  (using  Vqm(R)),  as  shown 
in  Figure  47.  Another  troubling  observation  is  that  some  of  the  probabilities  for  the 
lighter  elements  are  well  above  1,  which  is  of  course  mathematically  nonsense  (KHe 
using  VrgV)  is  the  worst,  with  a  probability  of  ~  500). 

Some  of  the  more  promising  observations  include  the  fact  that  all  of  the  meth¬ 
ods  of  calculation  that  were  presented  in  this  thesis  have  the  same  general  curve  of 
Ln(Probability)  vs.  Adiabaticity.  Another  positive  observation  is  that  calculations 
using  the  analytical  approximation  of  the  probability  integral  match  the  numerical 
calculations  with  Vrg(^)  very  well,  supplying  cross-validation  for  both  methods.  The 
analytic  approximation  offers  the  huge  benefit  of  being  able  to  make  quick  modifica¬ 
tions  to  the  radial  coupling  function’s  fit  parameters  in  Mathematica,  and  having  the 
results  show  up  almost  instantaneously.  Next,  the  calculations  using  VrgV)  match 
fairly  well  with  calculations  using  Vrg(R)-  The  approximations  made  in  converting 
Vrg{R)  to  VRcif)  appear  to  only  result  in  a  relatively  small  offset,  without  having 
any  significant  effect  on  the  shape  of  the  probability  curve.  The  use  of  the  quantum 
mechanical  coupling  functions  Vqm(R),  actually  introduced  a  much  larger  error  offset 
as  compared  to  the  rectangular-Gaussian  approximations.  This  may  have  been  due 
to  the  different  implementation  used  for  this  set  of  coupling  functions,  namely  the 
energy  conservation  requirement  which  was  necessary  to  overcome  the  difficulty  of 
not  having  data  for  R  <  3  [au].  Adding  energy  conservation  equates  with  remov- 
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Figure  47.  Log(Probability)  vs  Adiabaticity,  raw  results  for  all  atomic  species  (no  scaling),  along  with  available  experimental 
data. 


ing  an  approximation,  and  would  therefore  be  expected  to  result  in  a  more  accurate 
solution.  Getting  a  less  accurate  solution,  as  was  seen  when  comparing  the  results 
to  experiment  data,  seems  counterintuitive.  The  radial  coupling  functions,  Vqm(R), 
themselves  appear  to  be  validated,  at  least  for  potassium  and  rubidium,  through 
quantum  mechanical  calculations  of  temperature-dependent  cross  sections  which  were 
converted  to  probabilities  using  the  quantum-defect  cross  sections. 

An  important  thing  to  note  is  that  the  results  appear  to  be  very  sensitive  to  the 
shape  of  the  radial  coupling  function,  V.  An  area  of  further  research  could  be  to 
investigate  the  effect  of  different  trajectories,  which  would  effectively  be  changing 
the  shape  of  V(R(t)).  In  its  current  form,  the  C  code  developed  for  this  thesis  is 
able  to  perform  calculations  on  any  of  the  nine  M+Ng  pairs  presented  here  within 
ten  minutes  on  a  standard  desktop  computer,  with  calculations  on  some  pairs  taking 
only  seconds.  The  code  is  fairly  modular,  with  V (R)  and  R(b,  v,  t )  being  separate 
functions.  It  would  be  fairly  straightforward  to  introduce  different  forms  of  V  or  R 
into  this  code,  but  care  must  be  taken  if  the  form  of  V  used  is  no  longer  an  even 
function.  Approximations  were  made  assuming  that  V  was  even,  as  described  in 
Section  III. 


77 


Bibliography 


1.  Blank,  L,  D.  E.  Weeks,  and  G.  S.  Kedziora.  “M+Ng  potential  energy  curves 
including  spin-orbit  coupling  for  M  =  K,  Rb,  Cs  and  Ng  =  He,  Ne,  Ar”.  Journal 
of  Chemical  Physics ,  2012. 

2.  Curtis  F.  Gerald  and,  Patrick  O.  Wheatley.  Applied  Numerical  Analysis.  Pearson 
Education  Inc.,  7th  edition  edition,  2004. 

3.  Eshel,  Ben,  David  E.  Weeks,  and  Glen  P.  Perrarn.  “The  role  of  adiabaticity  in  alkali 
atom  fine  structure  mixing”.  SPIE  LASE ,  896207-896207.  International  Society 
for  Optics  and  Photonics,  2014. 

4.  Griffiths,  David  J.  Introduction  to  Quantum  Mechanics.  Pearson  Education,  Inc., 
2nd  edition,  2005. 

5.  Lewis,  Charlton  D.  Non-Adiabatic  Atomic  Transitions:  Computational  Cross  Sec¬ 
tion  Calculations  of  Alkali  Metal  -  Noble  Gas  Collisions.  Ph.D.  thesis,  Air  Force 
Institute  of  Technology,  September  2011. 

6.  Oppenheim,  Alan  V.  and  Ronald  W.  Schafer.  Discrete  Time  Signal  Processing. 
Pearson  Education  Inc.,  3rd  edition  edition,  2010. 

7.  Press,  William  H.,  Saul  A.  Teukolsky,  William  T.  Vetterling,  and  Brian  P.  Flan¬ 
nery.  Numerical  Recipes  in  C.  Cambridge  University  Press,  2nd  edition  edition, 
2002. 

8.  Weeks,  David  E.  “Fine  Structure  Transition  Probability  for  M+Ng  Collision  using 
the  Interaction  Function  Approximation  V{t)v .  Unpublished  Notes. 


78 


REPORT  DOCUMENTATION  PAGE 


Form  Approved 
OMB  No.  0704-0188 


Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std.  Z39.18 


